library(spatstat)
## Loading required package: spatstat.data
## Loading required package: spatstat.geom
## spatstat.geom 3.2-9
## Loading required package: spatstat.random
## spatstat.random 3.2-3
## Loading required package: spatstat.explore
## Loading required package: nlme
## spatstat.explore 3.2-7
## Loading required package: spatstat.model
## Loading required package: rpart
## spatstat.model 3.2-11
## Loading required package: spatstat.linnet
## spatstat.linnet 3.1-5
##
## spatstat 3.0-8
## For an introduction to spatstat, type 'beginner'
library(sp)
library(sf)
## Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE
library(splines)
load("wolf_clean_nondup.RData")
load("BC_Covariates.Rda")
#First moment descriptive statistics
plot(wolf_ppp, main='Coyotes')
## Quadracount
Q <- quadratcount(wolf_ppp, nx = 5, ny = 5)
plot(wolf_ppp,
pch = 20,
cex = 0.5,
main = "Quadrats Visualization")
plot(Q, add = TRUE, col = "red")
## Quadratest
quadrat.test(Q)
## Warning: Some expected counts are small; chi^2 approximation may be inaccurate
##
## Chi-squared test of CSR using quadrat counts
##
## data:
## X2 = 1319.3, df = 20, p-value < 2.2e-16
## alternative hypothesis: two.sided
##
## Quadrats: 21 tiles (irregular windows)
Q <- quadratcount(wolf_ppp, nx = 5, ny = 5)
plot(intensity(Q, image = TRUE), main = "Intensity", col = terrain.colors(256))
# Overlay the 'wolf_ppp' point pattern on the density plot
plot(wolf_ppp,
pch = 20, # Type of point, here '20' represents a filled circle.
cex = 0.6, # Point size.
col = "blue", # A single color for all points, you can choose your preferred color.
add = TRUE) # Indicates that points should be added to the existing plot.
## Heatmap
wolf_density <- density(wolf_ppp)
# Plot the density estimate
plot(wolf_density, main = "Heat map")
plot(wolf_ppp,
pch = 20, # Type of point, here '20' represents a filled circle.
cex = 0.5, # Point size.
col = "green", # Color of points, choose your preferred color.
add = TRUE) # Indicates that points should be added to the existing plot.
From above preliminary analysis, we find the distribution is
inhomogeneous.
Elevation <- BC_wolf$Elevation
Forest<-BC_wolf$Forest
HFI <- BC_wolf$HFI
Dist_Water <-BC_wolf$Dist_Water
We cut the 4 covariates into 5 classes and explore their relationship with coyote distribution.
elev_cut <- cut(Elevation,
5,
labels = c("very low", "low","medium","high","very high"))
# plot the elevation class image
plot(elev_cut,
main = "Elevation classes")
# overlay the park locations
plot(wolf_ppp,
cols = "black",
pch = 16,
cex = 0.6,
add = TRUE)
plot(wolf_ppp,
cols = "white",
pch = 16,
cex = 0.5,
add = TRUE)
wolf_elev_class <- elev_cut[wolf_ppp]
# Create a table of counts within each class
table(wolf_elev_class)
## wolf_elev_class
## very low low medium high very high
## 233 108 3 0 0
forest_cut <- cut(Forest,
5,
labels = c("very low", "low","medium","high","very high"))
# plot the elevation class image
plot(elev_cut,
main = "Forest classes")
# overlay
plot(wolf_ppp,
cols = "black",
pch = 16,
cex = 0.6,
add = TRUE)
plot(wolf_ppp,
cols = "white",
pch = 16,
cex = 0.5,
add = TRUE)
wolf_forest_class <- forest_cut[wolf_ppp]
# Create a table of counts within each class
table(wolf_forest_class)
## wolf_forest_class
## very low low medium high very high
## 165 83 22 52 19
HFI_cut <- cut(HFI,
5,
labels = c("very low", "low","medium","high","very high"))
# plot the elevation class image
plot(HFI_cut,
main = "HFI classes")
# overlay
plot(wolf_ppp,
cols = "black",
pch = 16,
cex = 0.6,
add = TRUE)
plot(wolf_ppp,
cols = "white",
pch = 16,
cex = 0.5,
add = TRUE)
wolf_HFI_class <- HFI_cut[wolf_ppp]
# Create a table of counts within each class
table(wolf_HFI_class)
## wolf_HFI_class
## very low low medium high very high
## 19 38 73 61 127
Dist_Water_cut <- cut(Dist_Water,
5,
labels = c("very close", "close","average","far","very far"))
# plot the elevation class image
plot(Dist_Water_cut,
main = "Dist_Water classes")
# overlay
plot(wolf_ppp,
cols = "black",
pch = 16,
cex = 0.6,
add = TRUE)
plot(wolf_ppp,
cols = "white",
pch = 16,
cex = 0.5,
add = TRUE)
Dist_Water_class <- Dist_Water_cut[wolf_ppp]
# Create a table of counts within each class
table(Dist_Water_class)
## Dist_Water_class
## very close close average far very far
## 332 13 2 0 0
rho_elev <- rhohat(wolf_ppp, Elevation)
rho_forest <- rhohat(wolf_ppp, Forest)
rho_hfi <- rhohat(wolf_ppp, HFI)
rho_water <- rhohat(wolf_ppp, Dist_Water)
plot(rho_elev, xlim=c(0, max(Elevation)), main="ρ vs Elevation")
Lower elevations have a higher than average density of coyote
occurrences.
plot(rho_forest, xlim=c(0, max(Forest)), main="ρ vs Forest")
When forest cover is <50%, the density seems higher.
plot(rho_hfi, xlim=c(0, max(HFI)), main="ρ vs HFI")
The intensity of sightings is initially low and increases sharply as the
HFI approaches 1.
plot(rho_water, xlim=c(0, max(Dist_Water)), main="ρ vs Dist_Water")
The relationship does not seem to be linear or simply exponential. The
significant width of the confidence interval in the latter part of the
plot for the distance from water indicates there is considerable
uncertainty in the intensity estimates at greater distances.
From above, it seems all covariates have different relationship with coyote distribution.
W <- owin(xrange = c(min(wolf_clean_unique$X), max(wolf_clean_unique$X)), yrange = c(min(wolf_clean_unique$Y), max(wolf_clean_unique$Y)))
wolf_ppp_2 <- ppp(x = wolf_clean_unique$X , y = wolf_clean_unique$Y, window = W)
miplot(wolf_ppp_2, main = "", pch = 16, col = "#046C9A")
Zoom in
miplot(wolf_ppp_2, main = "", xlim = c(0, 60000), pch = 16, col = "#046C9A")
Under the assumption of homogeneity, we can observe clear clustering
phenomena. However, this should be biased because inhomogeneity make
more sense.
#Estimate the empirical k-function
k_wolf <- Kest(wolf_ppp,correction=c("isotropic"))
k_wolf
## Function value object (class 'fv')
## for the function r -> K(r)
## ................................................................
## Math.label Description
## r r distance argument r
## theo K[pois](r) theoretical Poisson K(r)
## iso hat(K)[iso](r) Ripley isotropic correction estimate of K(r)
## ................................................................
## Default plot formula: .~r
## where "." stands for 'iso', 'theo'
## Recommended range of argument r: [0, 341660]
## Available range of argument r: [0, 341660]
plot(k_wolf,
main = "",
lwd = 2)
Bootstrapping with 95% CI
E_wolf <- envelope(wolf_ppp,
Kest,
correction="isotropic",
rank = 1,
nsim = 19,
fix.n = T)
## Generating 19 simulations of CSR with fixed number of points ...
## 1, 2, [4:30 remaining] 3,
## [4:04 remaining] 4, [3:47 remaining] 5, [3:31 remaining] 6,
## [3:21 remaining] 7, [3:04 remaining] 8, [2:49 remaining] 9,
## [2:33 remaining] 10, [2:19 remaining] 11, [2:03 remaining] 12,
## [1:48 remaining] 13, [1:33 remaining] 14, [1:17 remaining] 15,
## [1:01 remaining] 16, [46 sec remaining] 17, [30 sec remaining] 18,
## [15 sec remaining]
## 19.
##
## Done.
plot(E_wolf,
main = "",
lwd = 2)
For homogeneous proecess, it looks there is clustering.
lambda_wolf <- density(wolf_ppp, bw.ppl)
Kinhom_wolf <- Kinhom(wolf_ppp, lambda_wolf,correction="isotropic")
Kinhom_wolf
## Function value object (class 'fv')
## for the function r -> K[inhom](r)
## ................................................................................
## Math.label
## r r
## theo {K[inhom]^{pois}}(r)
## iso {hat(K)[inhom]^{iso}}(r)
## Description
## r distance argument r
## theo theoretical Poisson K[inhom](r)
## iso Ripley isotropic correction estimate of K[inhom](r)
## ................................................................................
## Default plot formula: .~r
## where "." stands for 'iso', 'theo'
## Recommended range of argument r: [0, 341660]
## Available range of argument r: [0, 341660]
# visualise the results
plot(Kinhom_wolf,
theo ~ r,
main = "",
col = "grey70",
lty = "dashed",
lwd = 2)
plot(Kinhom_wolf,
iso ~ r,
col = c("#046C9A"),
lwd = 2,
add = T)
Bootstrapping with 95% CI
lambda_wolf_pos <- density(wolf_ppp,
sigma=bw.ppl,
positive=TRUE)
#Simulation envelope (with points drawn from the estimated intensity)
E_wolf_inhom <- envelope(wolf_ppp,
Kinhom,
simulate = expression(rpoispp(lambda_wolf_pos)),
correction="border",
rank = 1,
nsim = 19)
## Generating 19 simulations by evaluating expression ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18,
## 19.
##
## Done.
# visualise the results
par(mfrow = c(1,2))
plot(E_wolf_inhom,
main = "",
lwd = 2)
# Zoom in on range where significant deviations appear
plot(E_wolf_inhom,
xlim = c(0,50000),
main = "",
lwd = 2)
After incorporating the assumption of inhomogeneity, significant
clustering only appears to exist within around 30000 meters.
pcf_wolf <- pcf(wolf_ppp)
pcf_wolf
## Function value object (class 'fv')
## for the function r -> g(r)
## ..............................................................
## Math.label Description
## r r distance argument r
## theo g[Pois](r) theoretical Poisson g(r)
## trans hat(g)[Trans](r) translation-corrected estimate of g(r)
## iso hat(g)[Ripley](r) isotropic-corrected estimate of g(r)
## ..............................................................
## Default plot formula: .~r
## where "." stands for 'iso', 'trans', 'theo'
## Recommended range of argument r: [0, 341660]
## Available range of argument r: [0, 341660]
plot(pcf_wolf)
# visualise the results
plot(pcf_wolf,
theo ~ r,
ylim = c(0,20),
main = "",
col = "grey70",
lwd = 2,
lty = "dashed")
plot(pcf_wolf,
iso ~ r,
col = c("#046C9A"),
lwd = 2,
add = T)
# Estimate g corrected for inhomogeneity
g_inhom <- pcfinhom(wolf_ppp)
# visualise the results
plot(g_inhom,
theo ~ r,
ylim = c(0,9),
main = "",
col = "grey70",
lwd = 2,
lty = "dashed")
plot(g_inhom,
iso ~ r,
col = c("#046C9A"),
lwd = 2,
add = T)
#Simulation envelope (with points drawn from the estimated intensity)
pcf_wolf_inhom <- envelope(wolf_ppp,
pcfinhom,
simulate = expression(rpoispp(lambda_wolf_pos)),
rank = 1,
nsim = 19)
## Generating 19 simulations by evaluating expression ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18,
## 19.
##
## Done.
# visualise the results
par(mfrow = c(1,2))
plot(pcf_wolf_inhom)
# Zoom in on range where significant deviations appear
plot(pcf_wolf_inhom,
xlim = c(0,60000),
main = "",
lwd = 2)
From the pair correlation function, it is evident that clustering is
more significant within 10,000 meters. Between 10,000 and 60,000 meters,
however, there is an indication of avoidance.
We check the collinearity of variables at first.
cor.im(BC_wolf$Elevation,BC_wolf$Forest,BC_wolf$HFI, BC_wolf$Dist_Water, use = "pairwise.complete.obs")
## ..1 ..2 ..3 ..4
## ..1 1.00000000 -0.26204718 -0.26625626 -0.03426387
## ..2 -0.26204718 1.00000000 0.06615541 0.04825439
## ..3 -0.26625626 0.06615541 1.00000000 0.13229408
## ..4 -0.03426387 0.04825439 0.13229408 1.00000000
From the correlation matrix above, we can see that the collinearity between Forest, HFI, Dist_Water and Elevation is very weak. All of them can be used for coyote spatial distribution modelling.
Based on the covariates relationships (rhohat plots), we guess that there is an non-linearity relationship between Dist_Water, HFI, Forest and Elevations. The quadratic term of Dist_water and forest might influence of spatial intensity of coyote, and the intensity of coyote might have strong correlation with HFI and Elevations’ higher polynomial terms.
We start at Forest, Dist_Water, and the Dist_Water and Elevation as our start point for ppm building. And utilize AIC as one of our comparsion standard
model_1 <- ppm(wolf_ppp ~ HFI + Forest + Dist_Water + Elevation, data = BC_wolf)
## Warning: Values of the covariate 'HFI' were NA or undefined at 0.67% (15 out of
## 2246) of the quadrature points. Occurred while executing: ppm.ppp(Q = wolf_ppp,
## trend = ~HFI + Forest + Dist_Water + Elevation,
model_1
## Nonstationary Poisson process
## Fitted to point pattern dataset 'wolf_ppp'
##
## Log intensity: ~HFI + Forest + Dist_Water + Elevation
##
## Fitted trend coefficients:
## (Intercept) HFI Forest Dist_Water Elevation
## -2.247316e+01 6.903940e+00 -2.058131e-03 9.045978e-06 -1.474970e-03
##
## Estimate S.E. CI95.lo CI95.hi Ztest
## (Intercept) -2.247316e+01 2.640138e-01 -2.299062e+01 -2.195570e+01 ***
## HFI 6.903940e+00 2.826821e-01 6.349893e+00 7.457986e+00 ***
## Forest -2.058131e-03 2.492556e-03 -6.943451e-03 2.827188e-03
## Dist_Water 9.045978e-06 2.622055e-05 -4.234535e-05 6.043731e-05
## Elevation -1.474970e-03 1.751079e-04 -1.818175e-03 -1.131765e-03 ***
## Zval
## (Intercept) -85.1211503
## HFI 24.4229792
## Forest -0.8257113
## Dist_Water 0.3449958
## Elevation -8.4232050
## Problem:
## Values of the covariate 'HFI' were NA or undefined at 0.67% (15 out of 2246)
## of the quadrature points
AIC(model_1)
## [1] 13741.68
The simple model shows that Forest and Dist_Water might have weak relationship with coyote’s distribution.
model_2 <- ppm(wolf_ppp ~ HFI + Forest + I(Forest^2) + I(Dist_Water) + Elevation, data = BC_wolf)
## Warning: Values of the covariate 'HFI' were NA or undefined at 0.67% (15 out of
## 2246) of the quadrature points. Occurred while executing: ppm.ppp(Q = wolf_ppp,
## trend = ~HFI + Forest + I(Forest^2) + I(Dist_Water) +
model_2
## Nonstationary Poisson process
## Fitted to point pattern dataset 'wolf_ppp'
##
## Log intensity: ~HFI + Forest + I(Forest^2) + I(Dist_Water) + Elevation
##
## Fitted trend coefficients:
## (Intercept) HFI Forest I(Forest^2) I(Dist_Water)
## -2.240181e+01 6.959465e+00 -1.012476e-02 9.729086e-05 -3.413309e-07
## Elevation
## -1.490738e-03
##
## Estimate S.E. CI95.lo CI95.hi Ztest
## (Intercept) -2.240181e+01 2.701828e-01 -2.293135e+01 -2.187226e+01 ***
## HFI 6.959465e+00 2.883950e-01 6.394222e+00 7.524709e+00 ***
## Forest -1.012476e-02 6.300284e-03 -2.247309e-02 2.223564e-03
## I(Forest^2) 9.729086e-05 6.939177e-05 -3.871451e-05 2.332962e-04
## I(Dist_Water) -3.413309e-07 2.736704e-05 -5.397975e-05 5.329709e-05
## Elevation -1.490738e-03 1.758718e-04 -1.835441e-03 -1.146036e-03 ***
## Zval
## (Intercept) -82.91350640
## HFI 24.13171099
## Forest -1.60703315
## I(Forest^2) 1.40205182
## I(Dist_Water) -0.01247233
## Elevation -8.47627596
## Problem:
## Values of the covariate 'HFI' were NA or undefined at 0.67% (15 out of 2246)
## of the quadrature points
AIC(model_2)
## [1] 13741.75
model_3 <- ppm(wolf_ppp ~ HFI + I(Dist_Water) + I(Forest^2) + I(Forest) + Elevation, data = BC_wolf)
## Warning: Values of the covariate 'HFI' were NA or undefined at 0.67% (15 out of
## 2246) of the quadrature points. Occurred while executing: ppm.ppp(Q = wolf_ppp,
## trend = ~HFI + I(Dist_Water) + I(Forest^2) +
model_3
## Nonstationary Poisson process
## Fitted to point pattern dataset 'wolf_ppp'
##
## Log intensity: ~HFI + I(Dist_Water) + I(Forest^2) + I(Forest) + Elevation
##
## Fitted trend coefficients:
## (Intercept) HFI I(Dist_Water) I(Forest^2) I(Forest)
## -2.240181e+01 6.959465e+00 -3.413309e-07 9.729086e-05 -1.012476e-02
## Elevation
## -1.490738e-03
##
## Estimate S.E. CI95.lo CI95.hi Ztest
## (Intercept) -2.240181e+01 2.701828e-01 -2.293135e+01 -2.187226e+01 ***
## HFI 6.959465e+00 2.883950e-01 6.394222e+00 7.524709e+00 ***
## I(Dist_Water) -3.413309e-07 2.736704e-05 -5.397975e-05 5.329709e-05
## I(Forest^2) 9.729086e-05 6.939177e-05 -3.871451e-05 2.332962e-04
## I(Forest) -1.012476e-02 6.300284e-03 -2.247309e-02 2.223564e-03
## Elevation -1.490738e-03 1.758718e-04 -1.835441e-03 -1.146036e-03 ***
## Zval
## (Intercept) -82.91350640
## HFI 24.13171099
## I(Dist_Water) -0.01247233
## I(Forest^2) 1.40205182
## I(Forest) -1.60703315
## Elevation -8.47627596
## Problem:
## Values of the covariate 'HFI' were NA or undefined at 0.67% (15 out of 2246)
## of the quadrature points
AIC(model_3)
## [1] 13741.75
model_4 <- ppm(wolf_ppp ~ HFI*I(Dist_Water) + Elevation, data = BC_wolf)
## Warning: Values of the covariate 'HFI' were NA or undefined at 0.67% (15 out of
## 2246) of the quadrature points. Occurred while executing: ppm.ppp(Q = wolf_ppp,
## trend = ~HFI * I(Dist_Water) + Elevation,
model_4
## Nonstationary Poisson process
## Fitted to point pattern dataset 'wolf_ppp'
##
## Log intensity: ~HFI * I(Dist_Water) + Elevation
##
## Fitted trend coefficients:
## (Intercept) HFI I(Dist_Water) Elevation
## -2.271192e+01 7.225534e+00 7.827459e-05 -1.492488e-03
## HFI:I(Dist_Water)
## -1.232840e-04
##
## Estimate S.E. CI95.lo CI95.hi Ztest
## (Intercept) -2.271192e+01 2.305704e-01 -2.316383e+01 -2.226002e+01 ***
## HFI 7.225534e+00 3.048343e-01 6.628070e+00 7.822998e+00 ***
## I(Dist_Water) 7.827459e-05 5.509886e-05 -2.971719e-05 1.862664e-04
## Elevation -1.492488e-03 1.764431e-04 -1.838310e-03 -1.146666e-03 ***
## HFI:I(Dist_Water) -1.232840e-04 9.712567e-05 -3.136468e-04 6.707884e-05
## Zval
## (Intercept) -98.503212
## HFI 23.703154
## I(Dist_Water) 1.420621
## Elevation -8.458754
## HFI:I(Dist_Water) -1.269324
## Problem:
## Values of the covariate 'HFI' were NA or undefined at 0.67% (15 out of 2246)
## of the quadrature points
AIC(model_4)
## [1] 13740.81
From those models above, we can see that all terms contains Dist_Water and Forest (no matter what kind it is) has week effect on coyote’s distribution. We remove them from our model.
model_5 <- ppm(wolf_ppp ~ HFI + Elevation, data = BC_wolf)
## Warning: Values of the covariate 'HFI' were NA or undefined at 0.67% (15 out of
## 2246) of the quadrature points. Occurred while executing: ppm.ppp(Q = wolf_ppp,
## trend = ~HFI + Elevation, data = list(list(
model_5
## Nonstationary Poisson process
## Fitted to point pattern dataset 'wolf_ppp'
##
## Log intensity: ~HFI + Elevation
##
## Fitted trend coefficients:
## (Intercept) HFI Elevation
## -22.579650812 7.012690936 -0.001477315
##
## Estimate S.E. CI95.lo CI95.hi Ztest
## (Intercept) -22.579650812 0.2064665368 -22.984317788 -22.174983836 ***
## HFI 7.012690936 0.2510511123 6.520639798 7.504742075 ***
## Elevation -0.001477315 0.0001747784 -0.001819874 -0.001134755 ***
## Zval
## (Intercept) -109.362278
## HFI 27.933320
## Elevation -8.452501
## Problem:
## Values of the covariate 'HFI' were NA or undefined at 0.67% (15 out of 2246)
## of the quadrature points
AIC(model_5)
## [1] 13738.63
Now all terms are effective in our model. And we want to explore if their polynomial fits the model better. Since the higher-degree polynomial will easily make the model fail to converge, we use gam instead of ppm
model_6 <- ppm(wolf_ppp ~ bs(HFI, knots = c(0.2,0.6), degree = 3, df = 4) + bs(Elevation, degree = 3, df = 3), data = BC_wolf, use.gam = TRUE)
## Warning: Values of the covariate 'HFI' were NA or undefined at 0.67% (15 out of
## 2246) of the quadrature points. Occurred while executing: ppm.ppp(Q = wolf_ppp,
## trend = ~bs(HFI, knots = c(0.2, 0.6), degree = 3,
model_6
## Nonstationary Poisson process
## Fitted to point pattern dataset 'wolf_ppp'
##
## Log intensity: ~bs(HFI, knots = c(0.2, 0.6), degree = 3, df = 4) +
## bs(Elevation, degree = 3, df = 3)
##
## Fitted trend coefficients:
## (Intercept)
## -23.8740545
## bs(HFI, knots = c(0.2, 0.6), degree = 3, df = 4)1
## 0.2301344
## bs(HFI, knots = c(0.2, 0.6), degree = 3, df = 4)2
## 5.1913730
## bs(HFI, knots = c(0.2, 0.6), degree = 3, df = 4)3
## 6.9626038
## bs(HFI, knots = c(0.2, 0.6), degree = 3, df = 4)4
## 6.5455994
## bs(HFI, knots = c(0.2, 0.6), degree = 3, df = 4)5
## 8.4684642
## bs(Elevation, degree = 3, df = 3)1
## -6.0917997
## bs(Elevation, degree = 3, df = 3)2
## 14.2016702
## bs(Elevation, degree = 3, df = 3)3
## -64.3914811
##
## For standard errors, type coef(summary(x))
## Problem:
## Values of the covariate 'HFI' were NA or undefined at 0.67% (15 out of 2246)
## of the quadrature points
AIC(model_6)
## [1] 13620.24
anova(model_5, model_6, test = "LRT")
## Warning: Models were re-fitted with use.gam=TRUE
## Warning: Values of the covariate 'HFI' were NA or undefined at 0.67% (15 out of
## 2246) of the quadrature points. Occurred while executing: ppm.ppp(Q = wolf_ppp,
## trend = ~HFI + Elevation, data = list(list(
## Warning: Values of the covariate 'HFI' were NA or undefined at 0.67% (15 out of
## 2246) of the quadrature points. Occurred while executing: ppm.ppp(Q = wolf_ppp,
## trend = ~bs(HFI, knots = c(0.2, 0.6), degree = 3,
## Analysis of Deviance Table
##
## Model 1: ~HFI + Elevation Poisson
## Model 2: ~bs(HFI, knots = c(0.2, 0.6), degree = 3, df = 4) + bs(Elevation, degree = 3, df = 3) Poisson
## Npar Df Deviance Pr(>Chi)
## 1 3
## 2 9 6 130.4 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
par_res_elev <- parres(model_6, "Elevation")
## Warning: Some infinite, NA or NaN increments were removed
par_res_hfi <- parres(model_6, "HFI")
## Warning: Some infinite, NA or NaN increments were removed
## Warning: Values for 1 query point lying outside the pixel image domain were
## estimated by projection to the nearest pixel
par_res_forest <- parres(model_6, "Forest")
## Warning: Some infinite, NA or NaN increments were removed
par_res_water <- parres(model_6, "Dist_Water")
## Warning: Some infinite, NA or NaN increments were removed
par(mfrow = c(1,2))
plot(par_res_elev,
legend = FALSE,
lwd = 2,
main = "",
xlab = "Elevation (m)")
plot(par_res_hfi,
legend = FALSE,
lwd = 2,
main = "",
xlab = "HFI")
As the AIC and anova report shown, our new GAM model has a better
performance now. But it also looks luck performance on the low
elevation, and the HFI’s confidence interval has high variance when HFI
near to 0.6 or lower than 0.2. We add two knots to capture more detailed
variations in coyote activity intensity with respect to elevation, and
we increased the weight proportion for HFI > 0.8 because, according
to rhohat, the frequency of coyote sightings only begins to rise sharply
after HFI exceeds 0.8.
model_7 <- ppm(wolf_ppp ~ bs(HFI, degree = 4, df = 4) + I(HFI > 0.8) + bs(Elevation, knots = c(300,1400),degree = 4, df = 3), data = BC_wolf, use.gam = TRUE)
## Warning: Values of the covariate 'HFI' were NA or undefined at 0.67% (15 out of
## 2246) of the quadrature points. Occurred while executing: ppm.ppp(Q = wolf_ppp,
## trend = ~bs(HFI, degree = 4, df = 4) +
model_7
## Nonstationary Poisson process
## Fitted to point pattern dataset 'wolf_ppp'
##
## Log intensity: ~bs(HFI, degree = 4, df = 4) + I(HFI > 0.8) + bs(Elevation,
## knots = c(300, 1400), degree = 4, df = 3)
##
## Fitted trend coefficients:
## (Intercept)
## -2.397222e+01
## bs(HFI, degree = 4, df = 4)1
## 2.502930e+00
## bs(HFI, degree = 4, df = 4)2
## 1.058737e+01
## bs(HFI, degree = 4, df = 4)3
## 3.644627e+00
## bs(HFI, degree = 4, df = 4)4
## 8.586450e+00
## I(HFI > 0.8)TRUE
## -5.246966e-02
## bs(Elevation, knots = c(300, 1400), degree = 4, df = 3)1
## 2.285253e-01
## bs(Elevation, knots = c(300, 1400), degree = 4, df = 3)2
## -4.159536e+00
## bs(Elevation, knots = c(300, 1400), degree = 4, df = 3)3
## 9.936512e+00
## bs(Elevation, knots = c(300, 1400), degree = 4, df = 3)4
## -3.594385e+01
## bs(Elevation, knots = c(300, 1400), degree = 4, df = 3)5
## 7.435627e+01
## bs(Elevation, knots = c(300, 1400), degree = 4, df = 3)6
## -1.106199e+03
##
## For standard errors, type coef(summary(x))
## Problem:
## Values of the covariate 'HFI' were NA or undefined at 0.67% (15 out of 2246)
## of the quadrature points
AIC(model_7)
## [1] 13604.38
anova(model_6, model_7, test = "LRT")
## Warning: Models were re-fitted with use.gam=TRUE
## Warning: Values of the covariate 'HFI' were NA or undefined at 0.67% (15 out of
## 2246) of the quadrature points. Occurred while executing: ppm.ppp(Q = wolf_ppp,
## trend = ~bs(HFI, knots = c(0.2, 0.6), degree = 3,
## Warning: Values of the covariate 'HFI' were NA or undefined at 0.67% (15 out of
## 2246) of the quadrature points. Occurred while executing: ppm.ppp(Q = wolf_ppp,
## trend = ~bs(HFI, degree = 4, df = 4) +
## Analysis of Deviance Table
##
## Model 1: ~bs(HFI, knots = c(0.2, 0.6), degree = 3, df = 4) + bs(Elevation, degree = 3, df = 3) Poisson
## Model 2: ~bs(HFI, degree = 4, df = 4) + I(HFI > 0.8) + bs(Elevation, knots = c(300, 1400), degree = 4, df = 3) Poisson
## Npar Df Deviance Pr(>Chi)
## 1 9
## 2 12 3 21.856 6.99e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
par_res_elev <- parres(model_7, "Elevation")
## Warning: Some infinite, NA or NaN increments were removed
par_res_hfi <- parres(model_7, "HFI")
## Warning: Some infinite, NA or NaN increments were removed
## Warning: Values for 1 query point lying outside the pixel image domain were
## estimated by projection to the nearest pixel
par_res_forest <- parres(model_7, "Forest")
## Warning: Some infinite, NA or NaN increments were removed
par_res_water <- parres(model_7, "Dist_Water")
## Warning: Some infinite, NA or NaN increments were removed
par(mfrow = c(1,2))
plot(par_res_elev,
legend = FALSE,
lwd = 2,
main = "",
xlab = "Elevation (m)")
plot(par_res_hfi,
legend = FALSE,
lwd = 2,
main = "",
xlab = "HFI")
Now We can see that the whole plot fitting has improved.
Since in the coyote intensity distribution, we can see that the coyote clustered in the HFI > 0.8. So we see the fitting in the small window.
par(mfrow = c(1,2))
plot(par_res_elev,
legend = FALSE,
lwd = 2,
main = "",
xlab = "Elevation (m)", ylim = c(-20,20))
plot(par_res_hfi,
legend = FALSE,
lwd = 2,
main = "",
xlab = "HFI", xlim = c(0.8,1))
We can see that the model fitting well now.
plot(model_7,log=TRUE,se = FALSE,superimpose = FALSE, n= 100)
Here is our residual plot
res <- residuals(model_7)
## Warning: Some infinite, NA or NaN increments were removed
na_indexes <- which(is.na(res$val))
# If you need to exclude NA values explicitly
if (length(na_indexes) > 0) {
res <- res[-na_indexes]
plot(res, cols='transparent')
} else {
plot(res, cols='transparent')
}
From the graph above, we can see that our model can fit the coyote distribution well. # lucking variable check
lurking(model_7, BC_wolf$Forest)
## Warning in LurkEngine(object = object, type = type, cumulative = cumulative, :
## 19 out of 2246 quadrature points discarded because they lie outside the domain
## of the covariate image
## Warning: model was fitted by gam(); asymptotic variance calculation ignores
## this
lurking(model_7, BC_wolf$Dist_Water)
## Warning: model was fitted by gam(); asymptotic variance calculation ignores
## this
There is less lucking variable effect in our model. So it means our
variable selected is suitable.
diagnose.ppm(model_7)
## Warning: Some infinite, NA or NaN increments were removed
## Warning: model was fitted by gam(); asymptotic variance calculation ignores
## this
## Warning: model was fitted by gam(); asymptotic variance calculation ignores
## this
## Model diagnostics (raw residuals)
## Diagnostics available:
## four-panel plot
## mark plot
## smoothed residual field
## x cumulative residuals
## y cumulative residuals
## sum of all residuals
## sum of raw residuals in entire window = -5.324e-10
## area of entire window = 9.483e+11
## quadrature area = 9.388e+11
## range of smoothed field = [-1.651e-10, 3.034e-10]
We have found an interesting issue that our residuals increase on the southwest coast, where we also have a higher record of coyote observations. Diagnostically, we can see that the peaks of the errors almost all occur in the urban areas of British Columbia: Vancouver, Kelowna, etc. We believe that these prediction errors stem from data sampling. Most of the coyote data relies on human sighting reports, which causes areas with high Human Footprint Index (HFI) to have more coyote records, while areas with less human activity have fewer records, hence our data exhibits sampling bias. Additionally, due to factors like urban environments and mountains, the distribution of coyotes is not homogeneous. This means that our four variables: forest, distance to water, elevation, and HFI, cannot fully encapsulate the factors influencing coyote distribution. To accurately predict coyote appearances, we need more data.